An efficient procedure for the development of optimized Projector Augmented Wave 

basis functions 
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In the Projector Augmented Wave (PAW) method, a local potential, basis functions, and pro- 
jector functions form an All-Electron (AE) basis for valence wave functions in the application of 
Density Functional Theory (DFT) . The construction of these potentials, basis functions and pro- 
jector functions for each element can be complex, and several codes capable of utilizing the PAW 
method have been otherwise prevented from its use by the lack of PAW basis sets for all atoms. We 
have developed a procedure that improves the ease and efficiency of construction of PAW basis sets. 
An evolutionary algorithm is used to optimize PAW basis sets to accurately reproduce scattering 
properties of the atom and which converge well with respect to the energy cutoff in a planewave 
basis. We demonstrate the procedure for the case of Ga with the 4s, 4p, and 3d electrons treated 
as valence. Calculations with this Ga PAW basis set are efficient and reproduce results of linearized 
augmented plane wave (LAPW) calculations. We also discuss the relationship between total energy 
convergence with respect to the energy cutoff and the magnitude of the matching radius of the PAW 
set. 



I. INTRODUCTION 



The Projector Augmented Wave (PAW) method of 
Blochl 1 provides an efficient, all-electron (AE) basis 
for Density Functional Theory 2 -^ (DFT) calculations 
within the frozen-core approximation. The method is 
applicable to all atoms, and the efficient treatment of 
first-row and transition metal elements compared with 
norm-conserving pseudopotential methods is improved 
dramatically^. The PAW method requires the construc- 
tion of a local potential, atom-centered radial (pseudo- 
and AE) basis functions, and projector functions for each 
atom. For PAW functions highly optimized with regard 
to computational efficiency, this construction can be com- 
plex, including a local potential and multiple sets of ba- 
sis and projector functions for each angular momentum 
channel and principal quantum number. To accelerate 
the construction of these PAW functions, and, further, 
to optimize with respect to the energy cutoff, we have 
developed an efficient, partially automated procedure for 
searching this high-dimensional parameter space. 

In section [TT| we introduce the basic elements of the 
PAW and its associated parameters. Our procedure for 
obtaining efficient ab initio PAW functions is developed 
in section IIII1 In section IIVI results of the construction 
of the Gallium (Ga) PAW are given including logarith- 
mic derivatives, total energy convergence, optimization 
results, and lattice constant and bulk modulus calcu- 
lations. A primary key to obtaining optimal total en- 
ergy convergence is a large matching radius, and this is 
shown by the comparison of a series of Ga PAW basis sets 
with increasing matching radii. A moderate variation of 
lattice constant with the matching radius parameter is 
found, and this will be discussed later. PAW basis sets 
constructed in this way are expected to accurately and ef- 
ficiently calculate the properties of solids within the DFT 
framework and the frozen-core approximation. 
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FIG. 1. (Color online) Ga AE and PAW local potentials. 

II. ELEMENTS OF THE PAW CONSTRUCTION 

The construction of PAW basis sets 5 - is similar to the 
construction of pseudopotentials. First, a reference con- 
figuration is chosen and electrons are divided into core 
and valence electrons. A local potential is created which 
is smooth within the PAW matching radius Rpaw and 
matches the AE potential at and beyond Rpaw, as in fig- 
ure!]] Radial basis functions and dual projector functions 
used as the basis for electronic wave functions are con- 
structed for each angular momentum channel such that 
a projection operation of p on a smooth (pseudo) basis 
function <p reproduces the AE basis function, <fi, as in fig- 
ure [3J A single basis function per angular momentum 
channel is often sufficient, but more regularly two func- 
tions per channel are required to approach a complete 
basis set^, and in rare cases multiple basis functions in 
one angular momentum channel but with different prin- 
cipal quantum numbers may also be needed to represent 
the wave functions of semi-core electrons. 

We use the 'atompaw' program 5 to construct basis sets 
in an RRKJ (Rappe, Rabe, Kaxiras, and Joannopou- 
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FIG. 2. (Color online) Ga Projector and Basis Functions. For 
each angular momentum channel, two basis functions may be 
used. The first PAW AE basis function, pseudo basis function, 
and projector function are shown in column a and the second 
set of functions is in column b. 



los) style 7 , one of several available in the 'atompaw' 
program. The RRKJ scheme optimizes the kinetic en- 
ergy because this dominates the total energy conver- 
gence. With RRKJ pseudo- basis functions and projec- 
tors, each projector function is associated with an indi- 
vidual matching radius, Rl,i, and an energy Elj. The 
energy of the first pseudo basis function in each angu- 
lar momentum channel is taken as the eigenvalue of the 
corresponding AE solution in the atomic system, leaving 
three free construction parameters per angular momen- 
tum channel when two projectors are used. 

A typical PAW will have 3 * Nl + 2 construction pa- 
rameters, including energy and matching radius param- 
eters for the local potential, where Nl is the number 
of angular momentum channels used in the basis, usu- 
ally corresponding to the angular momentum channels 
of the valence electrons. These parameters are chosen to 
optimize the matching of logarithmic derivatives and the 
total energy convergence with respect to an energy cutoff 
of the planewave basis when using these PAW functions 
to perform calculations of solids. Basis functions associ- 
ated with separate angular momentum channels appear 
to be independent with respect to logarithmic deriva- 
tives. However, in a crystalline calculation, these were 
found to have a strong interdependence in their effect on 
total energy convergence and in the self-consistent en- 
ergy minimization. Therefore, with a typical resolution 
of 0.05 Bohr in the matching radii and 0.1 Rydbergs in 
energy, and total ranges of 1 Bohr in the radii and 10 Ry- 
dbergs in energies, this creates an entire parameter space 
of 20 2Nl+1 * 100 Nl+1 total number of possible parame- 
ter sets, or about 1.28 * 10 17 sets for elements with s, p, 
and d basis functions. The size of this parameter space 
precludes a purely brute-force attempt to optimize PAW 
projection basis sets. In addition to previous advice in 
the literature and in various groups 8 " regarding the selec- 
tion of construction parameters, below and in section Hill 



we discuss further ways that this parameter space can be 
reduced and efficiently searched. 

The matching of logarithmic derivatives of the pseudo- 
and AE wave functions, evaluated with respect to the 
radius at Rpaw &s a function of energy, represents an 
equivalence of the scattering properties of the PAW basis 
set and the atom. The matching of logarithmic deriva- 
tives over a wide range of energies is one indication of a 
transferable basis set. In this ab initio construction of 
PAW basis sets, the primary aim is to accurately repro- 
duce logarithmic derivatives^, and we do not optimize 
the sets by the matching of results of calculations with 
any observable. 

Traditionally PAW parameters are chosen by hand, 
logarithmic derivatives are inspected visually, and post- 
testing of basis sets for optimal total energy convergence, 
and eventually for completeness by comparing with other 
AE results, is done individually for each constructed 
PAW. This may take many iterations to discover well- 
matching logarithmic derivatives which are also optimal 
for total energy convergence. Also, much of parame- 
ter space invariably remains unexplored due to its pro- 
hibitive size. In this paper we propose an approach to re- 
duce the time and labor associated with producing PAW 
basis sets and to discover otherwise unobtainable opti- 
mized PAWs. In our procedure, first logarithmic deriva- 
tives are given a visual inspection to define a limited 
range of possible construction parameters, thereby reduc- 
ing the available parameter space significantly. Then an 
evolutionary algorithm is used to further optimize both 
the matching of logarithmic derivatives and the total en- 
ergy convergence of a PAW. In some cases, it is nec- 
essary to also optimize with respect to the number of 
self-consistent iterations. 



III. METHOD 

The method described in this paper for obtaining an 
optimized PAW basis set consists of, first, the selection 
of configurational parameters, including the division of 
electrons into core and valence electrons and the selec- 
tion of the basic types of local potentials and basis func- 
tions. Next, parameter ranges are reduced through a 
visual inspection of logarithmic derivative matching and 
basis function smoothness. Then, an evolutionary al- 
gorithm further optimizes the PAW within the reduced 
paramter space. Lastly, a few tests are performed to en- 
sure the basic accuracy and tranferability of the PAW. 
Further testing of the PAW may be necessary depending 
on the intended calculation^. 

Initial parameter ranges can be restricted as follows. A 
local potential matching radius Hp aw should be as large 
as possible without resulting in sphere overlap in crys- 
talline, molecular, or molecular dynamics calculations 8 , 
in order to optimize the total energy convergence, as will 
be discussed in section HV.51 This is usually a little less 
than half the nearest neighbor distance in a particular 
system, with room to allow for ionic motion or relax- 
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ation. In our case, a Troullier-Martins style local po- 
tential is used, and the energy parameter often ranges 
from -2 to 5 Rydbergs, while energy parameters for the 
RRKJ projector functions have ranged from -2 to 8 Ry- 
dbergs. After the local potential matching radius Rpaw 
is determined, all individual matching radii Rl,% are set 
equal to this. Then a range for the local potential energy 
is determined by viewing the logarithmic derivative for 
the L = L max + 1 channel. With a suitable local po- 
tential energy chosen in this range, energies of the pro- 
jector functions are varied simultaneously to find regions 
with good logarithmic derivatives, but with a secondary 
goal of obtaining smooth basis functions as well. If well- 
matched logarithmic derivatives are unobtainable with 
only these energy parameters, either the Rpaw or the 
Rh,i may then also be altered to improve the logarith- 
mic derivatives. Ranges of parameters which yield well- 
matched logarithmic derivatives are then used within the 
optimization program. In the case of Ga, E s ranged from 
3.0 to 7.0 Rydbergs, E p ranged from 0.0 to 7.0 Rydbergs, 
Ed ranged from 0.0 to 3.5 Rydbergs, and Ei oc ranged 
from 0.0 to 5.0 Rydbergs, leaving an effective parameter 
space for Ga of about five million sets. 

In order to automate the process of optimizing log- 
arithmic derivatives and total energy convergence, we 
define fitness scores for both of these properties which 
will be used in an optimization algorithm. For this pur- 
pose a numerical comparison of logarithmic derivatives is 
evaluated for each potential PAW set, according to the 
summation, XXrae - Vi.PAw) 2 1 \exp(-ABS(dy/dx)), 
where in this expression yi t AB an d Hi, paw are the all- 
electron and pseudo- wave function logarithmic deriva- 
tives as functions of energy, respectively. The purpose of 
the denominator here is merely to attenuate the effect of 
the difference between the AE and PAW values near a 
divergence in the logarithmic derivatives, as can be seen 
for example in figure [3J Scores for each angular momen- 
tum channel are summed and typically normalized by an 
average ideal matching score. 

If logarithmic derivatives are not within a reasonable 
tolerance (of about 0.3), a penalty factor is given as the 
fitness score and returned to the optimization routine 
without further testing. Ghost state o 11 ' 12 may be de- 
tected by the presence of a divergence in the PAW log- 
arithmic derivative where there is none in the AE loga- 
rithmic derivative, and this case is given a large penalty 
factor. If logarithmic derivatives are free of ghost states 
and the matching of PAW and AE derivatives is within 
tolerance, crystalline tests at two successive energy cut- 
offs are performed, and their energy difference, normal- 
ized by an ideal difference (about 0.01 Ry), is used as 
a convergence score. Logarithmic derivative matching 
is emphasized with stepped weighting, using a factor of 
10 for the weight function until within an ideal toler- 
ance (about 0.03) at which point the weight is reduced 
to 1, equivalent to the weight of the total energy conver- 
gence scores. Care must be taken so that the test with 
large energy cutoff is well-converged, in order to prevent 
a plateau from forming in the total energy convergence 



figure as a function of energy cutoff. The Socorro- 13 DFT 
program with PAW functionality is used for these tests. 
The normalized logarithmic derivative and total energy 
convergence scores are then summed to produce a final 
fitness score. 

The Design Analysis Kit for Optimization and Teras- 
cale Applications^ (DAKOTA) offers convenient han- 
dling of input parameters, a parallel computation frame- 
work, and an interface to several optimization techniques 
and packages. Comprehensive testing of DAKOTA algo- 
rithms for searching PAW parameter space, or perhaps 
new algorithms, could reduce the time requirements for 
discovering optimal PAW functions, but we find that the 
evolutionary algorithm of the coliny package^ is suffi- 
cient. In the coliny _ea algorithm, after a random group 
of initial PAW parameter sets are evaluated, future pa- 
rameter sets are chosen as combinations of well-scoring 
sets from the previous iteration. In addition to the cross- 
ing of well-scoring parameter sets, random mutations are 
also applied to original sets and to crossed sets, allowing 
further exploration into diverse areas in parameter space. 
In our practice, we typically use the 'two_point' method 
of crossing with a 'cross_over' rate of 0.9, a mutation rate 
of 0.25 with the 'offset_cauchy' mutation method, and a 
population size of 100 and maximum number of itera- 
tions anywhere from 3,000 to 25,000, depending upon 
the number of parameters involved. In the 'two_point' 
method, optimal parent sets from the current generation 
are crossed by the division of parent sets into three re- 
gions, and the middle section is taken from one parent 
and the end sections from another parent. In the 'off- 
set_cauchy' method for mutation, a random number is 
generated according to a cauchy distribution with a mean 
of 0, and this is added to a selected parameter. Since the 
optimization with DAKOTA is done with nominal test- 
ing, further tests of several of the top scoring PAW basis 
sets are used to better gauge the total energy conver- 
gence and to verify the transferability of the PAW. As 
a function of total energy cutoff, we look at the lattice 
constant, bulk modulus, and total energy convergence 
properties to evaluate PAW basis set optimization. 



IV. RESULTS 

IV. 1. Logarithmic Derivative Matching 

The use of an evolutionary algorithm to optimize the 
numerical matching of logarithmic derivatives results in 
most cases with a nearly identical matching of loga- 
rithmic derivatives in the associated angular momentum 
channels, a feat which is possible but painstaking when 
optimizing by hand. The case of the Ga PAW is shown in 
figure[3l Additionally, in many cases the matching of log- 
arithmic derivatives can be obtained without the use of 
matching radii Rl.i as optimization parameters in each 
angular momentum channel; these can often be set to the 
local potential matching radius, Rpaw- The energy of 
expansion for a second projector function in each channel 
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FIG. 3. (Color online) The AE logarithmic 
the Ga atom and its PAW approximation 
wave function matching and the scattering properties of the 
atom. Excellent agreement in all angular momentum channels 
L=0,l,2,3 is shown. 
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FIG. 4. DAKOTA total score optimization. The best score of 
each group of 100 parameter sets is plotted on the y-axis. A 
lower score represents better logarithmic derivative matching 
and better total energy convergence. 



IV. 3. Variability of convergence and completeness 
of the Ga PAW 



then becomes the primary parameter. This reduces the 
parameter space significantly and eliminates the problem 
of kinks in projector functions, saving time both in the 
preliminary inspection of logarithmic derivatives for the 
restriction of the parameter ranges and in the DAKOTA 
optimization. 



Total energy convergence rates and physical properties 
of PAW sets with the top ten DAKOTA scores in this 
trial run appear to be nearly identical, as can be seen in 
figuresElEl and[3 This is due partly to a large measure of 
similarity in parameter values of the optimized PAW sets, 
but also due to the effectiveness of the PAW method and 
an insensitivity of calculated properties such as lattice 
constants to the exact parameters of the construction. 



IV. 2. DAKOTA Total Score Optimization 

The final score used for a fitness evaluation is a com- 
bination of scores from the logarithmic derivative match- 
ing and from the total energy convergence, as discussed 
above in section IIIII Figure @] shows the scoring of a 
typical DAKOTA run of around 5,000 parameter sets for 
the Ga PAW, using only E Lji with L=0,l,2 and E Loc 
as DAKOTA parameters while keeping all Rl,i equal to 
Rpaw ■ If an appropriate normalization for both the log- 
arithmic derivatives and the total energy convergence is 
used, the fitness score should converge to around 1.0. In 
this particular optimization, an ideal total energy conver- 
gence of 1 mRy was used for normalizing the total energy 
score, the difference between the 25 Rydberg and 100 Ry- 
dberg cutoff calculations. An ideal logarithmic derivative 
matching of 0.01 was used to normalize the logarithmic 
derivative score. In this case, the total energy conver- 
gence normalization was too ambitious, and the optimal 
fitness score appears to converge to around 6.15, but the 
order of magnitude is acceptable. The PAW with the 
best score, known internally as Ga_Rp2.1et001_5095 was 
used in figures [TJ [2l and [3j with all matching radii set to 
2.1 Bohr, E s found to be 4.287 Rydbergs, E p found to 
be 4.647 Rydbergs, E d found to be 0.797 Rydbergs, and 
Eioc found to be 2.248 Rydbergs. The standard devia- 
tion for the top ten scoring sets are 0.36, 1.79, 0.03, and 
0.09 for E s , E p , Ed, and Eioc, respectively. 
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FIG. 5. Total energy convergence rates for the top ten results 
from a particular run are indistinguishable. Matching radii 
Kla and Rpaw are all equal to 2.1 Bohr. A soft Nitrogen 
PAW is used here so that the total energy convergence is 
determined by the Ga PAW. 

In the evaluation of the lattice parameter for Gallium 
Nitride (GaN) in the zinc-blende structure, the Rpaw 
parameter for Ga had a noticeable effect on the lat- 
tice constant, as can be seen in figure [8j with increasing 
matching radii leading to smaller lattice constants, but 
with a total range of only 0.5% in the lattice constant. 
We have similarly generated a series of N PAW sets with 
increasing Rpaw parameters, but this had negligible ef- 
fect on the GaN lattice constant. 
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Gallium Lattice Constant Convergence 

Top 10 results from Rp2.1et001, (zb) GaN with hard N 
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FIG. 6. GaN zinc-blende lattice constants as functions of 
energy cutoff are indistinguishable. Matching radii Rl,i and 
Rpaw are all equal to 2.1 Bohr. A hard Nitrogen PAW is 
used here so any variation in physical properties should be 
attributed to Ga PAW differences. 



IV. 4. Transferability 

Transferability in the PAW method is effective due 
to its AE treatment of wave functions and potentials 
within the core region, limited only by the complete- 
ness of basis functions, and by the inherent limitations of 
the frozen-core approximation and the DFT. To demon- 
strate the transferability of the Ga PAW, we have cal- 
culated lattice constants and bulk moduli for zinc-blende 
GaN, zinc-blende GaAs, and zinc-blende GaP. The values 
calculated in the local density approximation (LDA) 16 , 
shown in Table U are within 0.2% error in the lattice 
constants in comparison with linearized augmented plane 
wave (LAPW) values. Bulk moduli are within 1% error, 
except for zinc-blende GaP, which differs from LAPW 17 
values by nearly 2%. 



Gallium Nitride Bulk Modulus 
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FIG. 7. GaN zinc-blende bulk moduli as functions of energy 
cutoff are virtually indistinguishable. Matching radii Rl,i and 
Rpaw are all equal to 2.1 Bohr. A hard Nitrogen PAW is used 
here so any variation in physical properties can be attributed 
to Ga PAW differences. 
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FIG. 8. The GaN zinc-blende lattice constant varies with the 
Rpaw matching radius. The N PAW used in these calcu- 
lations has a matching radius of Rpaw = 1-1 Bohr so that 
there is no overlap of spheres for any of the included points. 
A separate DAKOTA optimization was performed for each 
value of Rpaw- The maximum difference in lattice constant 
from Rpaw=1-7 to Rpaw=2.5 is less than 0.5%. 



IV. 5. Rpaw as key to total energy convergence 

While the El{ parameters do affect total energy con- 
vergence, the dominant bottleneck parameter is the lo- 
cal potential matching radius, Rpaw- For a series of 
PAW basis sets with increasing Rpaw, figure [9] shows 
a clear correspondence between the magnitude of the 
matching radius and the total energy convergence. With 
a large matching radius, pseudo-basis functions and the 
local potential can be made smooth and therefore may 
be expanded in a small number of plane waves. A small 
amount of overlap of matching radii of atoms in a solid 
calculation can have negligible effects, but in general the 
spheres should not overlap, limiting the size of the match- 
ing radius and the total energy convergence. 

GaN Total Energy Convergence 

as a function of Gallium Rpaw 
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FIG. 9. (Color online) Correspondence between total energy 
convergence and the matching radius Rpaw- 



V. CONCLUSION 

The procedure described in this paper automates the 
process of generating optimized PAW basis sets using an 
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TABLE I. PAW calculations of zinc-blende lattice constants and bulk moduli of GaN, GaAs, and GaP were performed in the 
LDA with energy cutoffs of 100 Ry (although convergence is expected as in figures [5] and 0) . Each of these used a single Ga 
PAW, Ga_Rp2.1et001_5095. 



Crystal 


PAW ID 


a [Bohr] 


B [MBar] 


LAPW a [Bohr] 


LAPW B [Bohr] 


a % error 


B % error 


zinc-blende GaN 


N.febsec5953 


8.439 


2.03 


8.447 


2.05 


-0.09 


-0.9 


zinc-blende GaAs 


Arsenic_d_1058 


10.619 


0.76 


10.620 


0.77 


-0.02 


-0.7 


zinc-blende GaP 


P.cea_22169 


10.226 


0.91 


10.242 


0.89 


-0.16 


1.89 



a LAPW data from 'atompaw' websitail. 



evolutionary algorithm to efficiently search a large pa- 
rameter space. In the example of the Ga PAW, the result- 
ing PAW matches all-electron scattering properties and 
is as efficient as the construction method and the partic- 
ular element permits. Efficiency and efficacy of the PAW 
was confirmed in a handful of crystalline environments. 
This method will be of assistance in ongoing efforts to 
produce efficient PAW sets for various DFT codes. 



VI. ACKNOWLEDGEMENTS 

We acknowledge many discussions with Natalie 
Holzwarth and Marc Torrent, as well as for helpful 
guidance given on the atompaw website^. We would 
like to thank Alan Tackett, Greg Walker, and Rachael 
Hansel for an introduction to the DAKOTA program. 
Brian Adams helped in a critical way with details of the 
DAKOTA program. Sandia is a multiprogram labora- 
tory operated by Sandia Corporation, a Lockheed Mar- 
tin Company, for the United States Department of En- 
ergy's National Nuclear Security Administration under 
Contract No. DE-AC04-94AL85000. 



1 P. E. Blochl, Phys. Rev. B 50, 17953 (1994). 

2 P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964). 

3 W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965). 

4 G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999). 

5 N. A. W. Holzwarth, A. R. Tackett, and G. E. Matthews, 
Computer Physics Communications 135, 329 (2001), 
ISSN 0010-4655. 

6 M. Y. Chou, Phys. Rev. B 45, 11465 (1992). 

7 A. M. Rappe, K. M. Rabe, E. Kaxiras, and J. D. 
Joannopoulos, Phys. Rev. B 41, 1227 (1990). 

8 M. Torrent and N. Holzwarth, A user's guide for atompaw 
code (2007), http://www.abinit.org/downloads/PAW/ 
AtomPAW2Abinit-Mamial-html/atompawUG .pdf 

9 V. Heine, in Solid State Physics, edited by F. Seitz and 
D. Turnbull (Academic Press, Inc, 1970), vol. 24. 

10 A. Kiejna, G. Kresse, J. Rogal, A. De Sarkar, K. Reuter, 
and M. Scheffler, Phys. Rev. B 73, 035404 (2006). 

11 X. Gonze, R. Stumpf, and M. Scheffler, Phys. Rev. B 44, 
8503 (1991). 



X. Gonze, P. Kackell, and M. Scheffler, Phys. Rev. B 41, 
12264 (1990). 

Socorro is developed at Sandia National Laboratories and 
available from http://dft.sandia.gov/Socorro/ 
M. Eldred, S. Brown, B. Adams, D. Dunlavy, D. Gay, 
L. Swiler, A. Giunta, W. Hart, J.-P. Watson, J. Eddy, 
et al., Dakota, a multilevel parallel object-oriented frame- 
work for design optimization, parameter estimation, uncer- 
tainty quantification, and sensitivity analysis: Version 4-0 
users manual (2006), sandia Technical Report SAND2006- 
6337, Version 4.2, Updated November 2008. 
W. Hart, The coliny optimization library (2004), 
http : //software . sandia. gov/Acro/Coliny 
J. P. Perdew and Y. Wang, Phys. Rev. B 45, 13244 (1992). 
N. Holzwarth et al., atompaw, 

http : //www. wfu. edu/~nat alie/papers/pwpaw/man. html 



